library(tidyverse)
library(ggpubr)
library(patchwork)

# flowjo output ------------
mc38_igg_14d <- read_csv('00_util_scripts/data/14d-post-mc38-Igg.txt.csv')

mc38_igg_14d <- mc38_igg_14d |>
  select(-4) |>
  set_names(c('name','MFI','prop')) |>
  slice_head(n = 33) |>
  mutate(group = str_extract(name, 'ii|it|tt|mock|ova|unstain'),
         group = case_match(group,
                            'ii' ~ 'II',
                            'it' ~ 'IT',
                            'tt' ~ 'TT',
                            'ova' ~ 'OVA',
                            .default = group),
         group = fct_relevel(group, 'OVA'),
         sex = str_extract(name, 'fem|mal'))

mc38_igg_14d |>
  filter((sex != 'fem' | group == 'OVA') & MFI < 3900) |>
  ggplot(aes(group, MFI, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 0.3) +
  geom_jitter(height = 0, width = 0.1) +
  stat_compare_means(comparisons = list(c('II','TT')), method = 't.test')+
  theme_pubr() +
  labs_pubr() +
  labs(y = 'MFI of IgG-binding MC38 cells')

g1 <- last_plot()

mc38_igg_14d |>
  filter((sex != 'fem' | group == 'OVA') & MFI < 3900) |>
  ggplot(aes(group, prop, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 0.3) +
  geom_jitter(height = 0, width = 0.1) +
  stat_compare_means(comparisons = list(c('II','TT')), method = 't.test')+
  theme_pubr() +
  labs_pubr() +
  labs(y = 'Percent of IgG-binding MC38 cells')

g2 <- last_plot()

g1 + g2 + patchwork::plot_annotation(title = 'Male mice at 14d post s.c. MC38 transplant')

## 36d tumor igg ---------
mc38_igg_36d <- read_csv('00_util_scripts/data/36d-mc38-igg.csv')

mc38_igg_36d <- mc38_igg_36d |>
  select(-4) |>
  set_names(c('name','MFI','prop')) |>
  slice_head(n = 32) |>
  mutate(group = str_extract(name, 'wt|it|tt|hel|Tube'),
         group = case_match(group,
                            'wt' ~ 'II',
                            'it' ~ 'IT',
                            'tt' ~ 'TT',
                            'hel' ~ 'OVA',
                            .default = group),
         group = fct_relevel(group, 'OVA'),
         sex = str_extract(name, 'female|_male'))

mc38_igg_36d |>
  filter((sex != '_male' | (group == 'OVA' & prop < 40))) |>
  ggplot(aes(group, MFI, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 0.3) +
  geom_jitter(height = 0, width = 0.1) +
  stat_compare_means(comparisons = list(c('II','TT')), method = 't.test')+
  theme_pubr() +
  labs_pubr() +
  labs(y = 'MFI of IgG-binding MC38 cells')

g1 <- last_plot()

mc38_igg_36d |>
  filter(sex != '_male' | (group == 'OVA' & prop < 40)) |>
  ggplot(aes(group, prop, color = group)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 0.3) +
  geom_jitter(height = 0, width = 0.1) +
  stat_compare_means(comparisons = list(c('II','TT')), method = 't.test')+
  theme_pubr() +
  labs_pubr() +
  labs(y = 'Percent of IgG-binding MC38 cells')

g2 <- last_plot()

g1 + g2 + plot_annotation(title = 'Female mice at 36d post s.c. MC38 transplant')

# base R string exercise -------
parse_phone_number <- function(number_string) {
  res <- gsub(x = number_string, pattern = '[^0-9]', replacement = '')
  
  if (nchar(res) == 11 & substr(res, 1, 1) == '1') {
    res <- sub(x = res, pattern = '^1', replacement = '')
  }
  
  if (nchar(res) == 10 & grepl(x = substr(res, 1, 1), pattern = '[2-9]') & grepl(x = substr(res, 4, 4), pattern = '[2-9]') ) {
    return(res)}
  
  else NULL
}

parse_phone_number('223.456.7890')

# pie chart --------
pie_chart <- tibble(lambda_2 = rpois(n = 2000, lambda = 2),
                    lambda_3 = rpois(n = 2000, lambda = 3),
                    lambda_4 = rpois(n = 2000, lambda = 4),
                    lambda_5 = rpois(n = 2000, lambda = 5)) |>
  pivot_longer(everything(),
               names_to = 'lambda',
               values_to = 'values')

pie_chart |>
  mutate(bins = ntile(values, n = 6)) |>
  group_by(bins) |>
  mutate(bins_range = range(values) |> str_flatten('-')) |>
  ungroup() |>
  group_by(lambda) |>
  count(bins_range) |>
  ggplot(aes(x = '', y = n, fill = bins_range)) +
  geom_col() +
  coord_polar(theta = 'y') +
  scale_fill_viridis_d() +
  facet_wrap(~ lambda)

# tumor size plot ----
readxl::read_excel('~/Documents/Lab/rs1050501/mice/2304-MC38-3e5.xlsx',
                   skip = 1) ->
  mc38

mc38 <- read_delim('00_util_scripts/data/230324-mc38-I232T.txt')

mc38 <- mc38 |>
  pivot_longer(4:last_col(),
               names_to = 'days',
               names_transform = list(days = as.numeric),
               values_to = 'volume',
               values_transform = list(volume = as.numeric)) |>
  mutate(days = case_match(days, 11 ~ 13, 13 ~ 11, .default = days)) 

mc38_vol_test <- mc38 |>
  select(sex, days, genotype, volume) |>
  pivot_wider(names_from = genotype, 
              values_from = volume,
              values_fn = list) |>
    rowwise() |>
    mutate(II_IT = t.test(unlist(II), unlist(IT))$p.value,
           II_TT = t.test(unlist(II), unlist(TT))$p.value,
           IT_TT = t.test(unlist(TT), unlist(IT))$p.value) |>
  compare_means(volume ~ genotype, ref.group = 'TT', method = 't.test')
  
mc38 |>
  filter(sex == 'male') |>
  ggplot(aes(days, volume, color = genotype)) +
  #geom_jitter(width = 0.3, height = 0) +
  stat_summary(geom = 'point', fun = 'mean', size = 3) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 1, linewidth = 1) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 2) +
  theme_pubr() +
  labs_pubr() +
  scale_color_manual(values = c('grey','blue','red'),
                     label = c('II (n=5)','IT (n=3)','TT (n=5)')) +
  facet_wrap(~sex, scales = 'free', ncol = 1) +
  ylab('tumor volume (mm^3)')

# tumor weight ----
mc38_weight <- read_delim('00_util_scripts/data/230324-mc38-final-weight.txt', col_types = c('cccd'))

mc38_weight |>
  ggplot(aes(genotype, weight, color = genotype)) +
  stat_summary(geom = 'col', fun = 'mean', fill = 'white') +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = 0.3) +
  geom_jitter(height = 0, width = 0.1) +
  stat_compare_means(comparisons = list(c('II','TT'), c('IT','II')),
                     method = 't.test')+
  theme_pubr() +
  labs_pubr() +
  labs(y = 'Weight of tumor (g)')+
  facet_wrap(~sex)

# mIHC of CRC ---------
mihc <- readxl::read_excel('CRC-I/data/mIHC分析结果.xlsx')

mihc |>
  dplyr::count(`FCGR2B-I232T`)

mihc <- mihc |>
  fill(`FCGR2B-I232T`)

mihc <- mihc |>
  set_names(c('FCGR2B_IT','id','location','CD68','PANCK','CD163+CD68','HLADR+CD68+'))

mihc |>
  filter(location == 'Stroma') |>
  ggplot(aes(FCGR2B_IT, CD68, color = FCGR2B_IT)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'CD68+ cells (MF) in tumor region mIHC',
       y = 'Fraction of cells') +
  stat_compare_means(comparisons = list(c('TT','II'),
                                        c('TT','IT')),
                     method = 't.test')

mihc |>
  filter(location == 'Stroma') |>
  ggplot(aes(FCGR2B_IT, PANCK, color = FCGR2B_IT)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'PANCK+ cells (tumor cell marker) in tumor region mIHC') +
  stat_compare_means(comparisons = list(c('TT','II'),
                                        c('TT','IT')),
                     method = 't.test')

mihc |>
  filter(location == 'Stroma') |>
  ggplot(aes(FCGR2B_IT, `HLADR+CD68+`, color = FCGR2B_IT)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'HLA-DR+CD68+ M1 macrophage in stroma region mIHC') +
  stat_compare_means(comparisons = list(c('TT','II'),
                                        c('TT','IT')),
                     method = 't.test')

mihc |>
  filter(location == 'Stroma') |>
  ggplot(aes(FCGR2B_IT, `CD163+CD68`, color = FCGR2B_IT)) +
  stat_summary(fun = 'mean', geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'CD163+CD68+ M2 macrophage in stroma region mIHC') +
  stat_compare_means(comparisons = list(c('TT','II'),
                                        c('TT','IT')),
                     method = 't.test')

# 5d protocol of AOM DSS -----
dss5 <- read_delim('00_util_scripts/data/aom-dss-5d-221218.txt') |>
  mutate(across(3:last_col(), \(x)x/`0`)*100) |>
  pivot_longer(3:last_col(),names_to = 'day',values_to = 'weight') |>
  mutate(day = str_remove(day,'day') |> as.numeric())

dss5 |>
  ggplot(aes(day, weight, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(y = 'Weight change (%)', color = 'Fcgr2b-I232T genotype') +
  scale_color_discrete(label = c('TT (n=9)','II (n=12)')) +
  theme_pubr() +
  labs_pubr()

g1 <- last_plot()

dss5b2 <- read_delim('00_util_scripts/data/aom-dss-5d-221218b2.txt') |>
  mutate(across(3:last_col(), \(x)x/`0`)*100) |>
  pivot_longer(3:last_col(),names_to = 'day',values_to = 'weight') |>
  mutate(day = str_remove(day,'day') |> as.numeric())

dss5b2 |>
  ggplot(aes(day, weight, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(y = 'Weight change (%)', color = 'Fcgr2b-I232T genotype') +
  scale_color_discrete(label = c('TT (n=9)','II (n=12)')) +
  theme_pubr() +
  labs_pubr()

dss5 |>
  bind_rows(dss5b2) |>
  ggplot(aes(day, weight, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(x = 'Day', y = 'Weight change (%)',
       color = 'Fcgr2b-I232T genotype') +
  scale_color_discrete(label = c('TT (n=9)','II (n=12)')) +
  theme_pubr() +
  labs_pubr()

dss5 |>
  bind_rows(dss5b2) |>
  filter(day > 0) |>
  pivot_wider(id_cols = -id, names_from = genotype, values_from = weight,
              values_fn = list) |>
  rowwise() |>
  mutate(pval = t.test(WT, TT)$p.value) |>
  arrange(pval)

## disease activity index -----------
dss5bs <- read_delim('00_util_scripts/data/aom-dss-5d-stool221218.txt') |>
  pivot_longer(3:last_col(),names_to = 'day',values_to = 'stool') |>
  mutate(day = str_remove(day,'day') |> as.numeric())

dss5_wide <- dss5 |>
  pivot_wider(names_from = day, values_from = weight)
  
dss5_dai <- 4:14 |>
  map(\(x)dss5_wide[,x] / dss5_wide[,x-1]) |>
  list_cbind() |>
  bind_cols(dss5_wide[,1:2]) |>
  pivot_longer(1:11,names_to = 'day',names_transform = as.numeric, values_to = 'weightloss') |>
  right_join(dss5bs) |>
  mutate(weightloss_score = case_when(
    weightloss < .82 ~ 4,
    weightloss < .9 ~ 3,
    weightloss < .95 ~ 2,
    weightloss < 1 ~ 1,
    .default = 0),
    dai = weightloss_score + stool
  )

dss5_dai |>
  ggplot(aes(day, dai, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(y = 'Disease activity index', color = 'Fcgr2b-I232T genotype') +
  theme_pubr() +
  labs_pubr()

dss5bsb2 <- read_delim('00_util_scripts/data/aom-dss-5d-stool221218b2.txt') |>
  pivot_longer(3:last_col(),names_to = 'day',values_to = 'stool') |>
  mutate(day = str_remove(day,'day') |> as.numeric())

dss5b2_wide <- dss5b2 |>
  pivot_wider(names_from = day, values_from = weight)

dss5b2_dai <- 4:14 |>
  map(\(x)dss5b2_wide[,x] / dss5b2_wide[,x-1]) |>
  list_cbind() |>
  bind_cols(dss5b2_wide[,1:2]) |>
  pivot_longer(1:11,names_to = 'day',names_transform = as.numeric, values_to = 'weightloss') |>
  right_join(dss5bsb2) |>
  mutate(weightloss_score = case_when(
    weightloss < .82 ~ 4,
    weightloss < .9 ~ 3,
    weightloss < .95 ~ 2,
    weightloss < 1 ~ 1,
    .default = 0),
    dai = weightloss_score + stool
  )

dss5b2_dai |>
  ggplot(aes(day, dai, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(y = 'Disease activity index', color = 'Fcgr2b-I232T genotype') +
  theme_pubr() +
  labs_pubr()

dss5b2_dai |>
  bind_rows(dss5_dai) |>
  ggplot(aes(day, dai, group = genotype, color = genotype)) +
  stat_summary(geom = 'path', fun = 'mean', linewidth = 1) +
  stat_summary(geom = 'errorbar', fun.data = 'mean_se', width = .4) +
  labs(x = 'Day', y = 'Disease activity index',
       color = 'Fcgr2b-I232T genotype') +
  theme_pubr() +
  labs_pubr()

g2 <- last_plot()

g1 / g2

dss5b2_dai |>
  bind_rows(dss5_dai) |>
  slice_max(day) |>
  pivot_wider(id_cols = day, names_from = genotype, values_from = dai,
              values_fn = list) |>
  rowwise() |>
  mutate(pval = t.test(WT, TT)$p.value)

## survival ----------
library(survival)
library(survminer)
dss_srv <- read_delim('00_util_scripts/data/aom-dss-5d-survival.txt')

survfit(Surv(day, dead) ~ genotype, data = dss_srv) |>
  ggsurvplot(data = dss_srv,
             xlab = "Days post AOM treatment",
             pval = TRUE,
             risk.table = F,
             break.x.by = 50,
             legend.labs = c('TT (n=9)','WT (n=12)'))

## sacrifice stat ------
cac.sacri <- read_delim('00_util_scripts/data/240616-cac-tumor.txt')

cac.sacri |>
  ggplot(aes(genotype, length, fill = genotype)) +
  stat_summary(geom = 'col') +
  stat_summary(geom = 'errorbar', width = .5) +
  geom_jitter(width = .3, height = 0) +
  stat_compare_means(method = 't.test') +
  theme_pubr()

cac.sacri |>
  mutate(number = ifelse(genotype == 'WT', number + 1, number)) |>
  ggplot(aes(genotype, number, fill = genotype)) +
  stat_summary(geom = 'col') +
  stat_summary(geom = 'errorbar', width = .5) +
  geom_jitter(width = .3, height = 0) +
  stat_compare_means(method = 't.test') +
  theme_pubr()

# apc-cko weight ----------
apc.cko <- read_delim('CRC-I/data/mm_l012/vil.cre.apc.flx.1224.tsv') |>
  mutate(date = ymd(date))

apc.cko.3wpt <- apc.cko |>
  pivot_longer(-1) |>
  mutate(group = ifelse(name == 'II', 'II', 'TT')) |>
  mutate(ref = min(value), .by = name) |>
  mutate(ratio = value/ref) |>
  slice_max(date) 

apc.cko.3wpt |>
  ggplot(aes(group, ratio*100, color = group)) +
  geom_boxplot(outliers = F) +
  geom_jitter(width = .05, height = 0) +
  expand_limits(y = 100) +
  theme_bw() +
  labs_pubr() +
  labs(title = 'Weight change of Vil-CreER Apc-flox mice 3 week post induction',
       y = 'Weight proportion to original (%)') +
  annotate(geom = 'segment', x = 1, xend = 2, y = 111) +
  annotate(geom = 'text', x = 1.5, y = 111.5, label = 'p=0.027')

t.test(c(.01,.04,.04,.05)*10)

# ARRB1-T370M T2D PheWAS -------------
arrb1_phewas <- read_csv('mission/PheWAS_associations.csv', skip = 1,
                         col_names = c('pos','abbr','pheno','group','dicto','chrom','ppos','p','beta','stderr','n'))

arrb1_phewas
arrb1_tidy <- arrb1_phewas |>
  mutate(pheno = str_remove(pheno, 'description:') |> fct_reorder(p, .desc = T),
         stderr = stderr * 1.96,
         p = str_c('p=',signif(p,2)),
         n = str_c('n=',n)) |>
  select(pheno, p, beta, stderr, n)

arrb1_tidy |>
  ggplot(aes(pheno, beta, ymax = beta + stderr, ymin = beta - stderr)) +
  geom_pointrange(color = 'red') +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  geom_text(aes(label = p, y = -.8), size = 5) +
  geom_text(aes(label = n, y = -.3), size = 5) +
  expand_limits(y = -1) +
  coord_flip() +
  theme_pubr() +
  labs_pubr() +
  theme(plot.title.position = 'plot') +
  labs(x = 'phenotype', y = 'Log odds ratio (Beta)',
       title = 'PheWAS significant results of ARRB1-T370M on T2D-associated traits')
  
# IVIS L-012 imaging -------
'DSS	3.94E+04
DSS	6.63E+04
DSS	4.85E+04
Ctrl	2.82E+04
Ctrl	1.58E+04' |> read_delim(col_names = c('group','count')) |>
  ggplot(aes(group, count, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(width = .01, height = 0) +
  stat_summary(fun.data = 'mean_se', geom = 'errorbar', width = .2) +
  theme_pubr() +
  scale_color_manual(values = c('blue', 'red')) +
  labs(y = 'Photon count') +
  stat_compare_means(comparisons = list(c('Ctrl', 'DSS')), method = 't.test')

ivis0419 <- read_delim('CRC-I/data/mm_l012/total-radiance.txt',name_repair = 'universal') 

ivis0419 |>
  mutate(group = if_else(str_ends(Image.Number, '27|26'), 'II', 'TT')) |>
  ggplot(aes(group, Total.Flux..p.s., color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(width = .01, height = 0) +
  stat_summary(fun.data = 'mean_se', geom = 'errorbar', width = .2) +
  theme_pubr() +
  scale_color_manual(values = c('blue', 'red'),
                     labels = c('II (n=11)','TT (n=8)')) +
  labs(y = 'Normalized radiance') +
  stat_compare_means(comparisons = list(c('II', 'TT')),
                     method = 't.test', color = 'black')

# exercism stat ------
"
Python 348k
JavaScript 289k
Java 138k
C# 87k
C++ 90k
C 74k
Go 94k
R 10k
Julia 15k
Bash 47k
" |> read_delim(col_names = c('language', 'user')) |>
  mutate(user = str_remove(user, 'k') |> as.double(),
         language = fct_reorder(language, user)) |>
  ggplot(aes(language, user, fill = language)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'none') +
  labs(title = 'Users of different language tracks on Exercism',
       y = 'User count (k)')

# mean pLDDT of AF3 -------
library(jsonlite)
f0 <- read_json('fold_human_trpm7_monomer_full_data_3.json', simplifyVector = T)
f0$atom_plddts |> mean()

